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MLTE 

AN  LTE  SOLVER  FOR  DIATOMIC  MOLECULES 
AND  ATOMS  ABOVE  6,000  °K 


1.  Introduction 


Consider  an  isolated  system  containing  a  mixture  of  various  species.  Under  the 
condition  of  thermodynamic  equilibrium  the  particle  number  of  each  species  N, ,  is  that  set 
of  numbers  which  is  most  probable  to  be  observed  in  the  system  subject  to  certain 
constraints.  These  constraints  are  conservation  of  two  of  the  following  three:  mass, 
volume,  or  pressure.  Thermodynamic  equilibria  is  appropriate  to  systems  at  high  pressure 
and  for  timescales  long  compared  to  l/(reaction  rate).  Such  conditions  generally  hold  in 
automobile  combustion  engines,  waste  treatment  plants,  and  atmospheric  arc  torches. 

Even  in  systems  which  are  not  in  thermodynamic  equilibrium,  such  as  low  pressure 
inductively  coupled  plasma  processing  reactors,  the  rapid  calculation  of  species 
composition  assuming  equilibria  provides  a  first  order  guess  to  further  kinetic  calculations. 
Furthermore,  one  often  finds  that  while  an  entire  collection  of  species  is  not  in 
thermodynamic  equilibria,  certain  subsets  of  the  species  are,  due  to  the  rapid  reaction  rates 
among  the  members  of  the  set.  This  state  is  often  referred  to  as  one  in  partial  equilibria. 

According  to  Besmann  [Ref.l],  there  have  been  two  basic  approaches  to  solving 
for  the  most  probable  set  of  abundances  [  N,).  Brinkley  [Ref. 2,3]  in  the  1940's 
determined  chemical  equilibrium  constants  for  each  chemical  equation,  while  White,  et  al., 
[Ref.4]  developed  in  1958  an  approach  using  Lagrange  multipliers  to  directly  minimize  the 
free  energy  of  the  system  and  obtain  thermodynamic  equilibria.  A  good  description  of 
this  technique  from  1968  is  presented  by  Zeleznik  and  Gordon  pief.5].  Since  this  time 
several  generations  of  "standard"  equilibria  solvers  for  use  on  mainframe  computers  have 
emerged,  among  them  are  SOLGASMIX  [Ref.l]  and  CHEMKIN  [Ref. 6].  To  evaluate 
the  free  energy  of  a  system  one  needs  the  thermodynamic  data  for  each  species  in  the 
mixture.  A  standard  reference  for  the  required  functions,  such  as  specific  heat  capacity,  is 
McBride,  Gordon,  and  Reno  [Ref. 7].  Today  there  are  readily  available  software  packages 
for  Personal  Computers  which  come  with  their  own  data  base  for  several  thousand 
compounds. 

A  study  of  the  tables  in  Ref.7  for  the  thermodynamic  quantities  used  in  minimizing 
the  free  energy  reveals  that  the  values  for  each  compound  listed  in  the  table  are  limited  to 
temperatures  less  than  6000  K.  This  upper  limit  is  certainly  reasonable  for  any  research  in 
combustion  processes.  However,  there  are  conditions  in  reactors  used  for  plasma 
processing  which  easily  exceed  this  limit.  For  instance,  Mehmetoglu  and  Gauvin  [Ref.8] 
report  peak  temperatures  of  18,000  °K  in  a  transferred  arc  torch  and  Chen,  et  al.,  [Ref. 9] 
have  measured  temperatures  of  14,000  °K  in  a  plasma  jet  (non-transferred)  using  absolute 
emissivities.  The  purpose  of  the  present  research  is  to  describe  a  numerical  code 
developed  in  Code  6720  which  can  solve  for  thermodynamic  equilibria  among  a  mixture  of 
elemental  atoms,  diatomic  molecules,  associated  ions,  and  electrons  up  to  arbitrary 
temperatures.  Of  course  at  temperatures  where  all  of  the  molecules  have  been  dissociated 

Manuscript  approved  May  1 1,  1998. 


1 


(-20,000  °  K)  and  all  that  remains  are  atoms,  ions,  and  electrons,  then  the  standard  Salia 
equilibrium  formula  is  valid.  Below  temperatures  of  complete  dissociation  it  is  also 
possible  to  include  molecules  within  a  StUia-like  framework  and  solve  a  set  of  Saha 
equations  which  involve  the  partition  function  for  each  species  considered  [Ref.lOJ.  In 
this  approach  to  obtaining  equilibrium,  the  proper  set  of  abundances  { A', )  is  the  one  which 
maximizes  the  system  partition  function. 

In  Section  II  the  theoretical  approach  used  in  the  equilibrium  solver  is  developed 
with  an  example  consisting  of  two  elements;  carbon  and  oxygen.  Tlie  numerical  code  can 
treat  an  arbitrary  number  of  elements,  but  again  we  note  that  it  is  limited  to  only 
diatomic  molecules  and  their  associated  ions.  Section  III  is  a  user's  guide  to  MLTE,  the 
numerical  code,  and  presents  results  for  four  examples.  MLTE  is  a  subroutine  and  a 
listing  of  a  sample  driver,  written  in  FORTRAN,  is  included  in  Section  IV.  One  major 
component  of  the  code  is  the  calculation  of  the  partition  function  of  each  species  at  an 
arbitrary  temperature  from  a  data  base.  The  entries  for  the  molecular  and  atomic  data 
bases  are  described  in  Section  V.  Basically  the  molecular  data  for  the  energy  levels  in 
diatomic  molecules  is  taken  from  Herzberg’s  classic  book  [Ref.l  1]  or  Huber  and  Herzberg 
[Ref.  12].  Based  on  the  data  tables  it  is  just  as  easy  to  calculate  a  specific  heat  as  a  partition 
function.  If  the  specific  heat  for  an  individual  species  can  be  fit  above  6()(X)K  by  a  quartic 
in  temperature  and  the  heat  of  formation  is  known,  then  the  coefficients  can  be  formulated 
in  a  manner  similar  to  that  in  the  NASA  tables  [Ref.7].  It  is  then  possible  to  use  standard 
solvers  such  as  CHEMKIN  to  properly  calculate  the  equilibrium  abundances  at 
temperatures  above  the  existing  table  limits  of  6000K.  This  approach  is  described  in 
Section  VI  with  several  other  examples. 
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11.  MLTE 

A  Molecular  (and  Atomic)  LTE  Populations  Code 

The  purpose  of  the  MLTE  code  is  to  calculate  the  LTE  population  distribution 
over  the  set  of  allowed  states,  both  molecular  (diatomic  only)  and  atomic.  The  sets  of 
electronic  states  are  defined  in  the  input  data  files.  The  code  calculates  their  LTE 
distribution  at  temperature  .  For  the  molecules  the  distribution  over  the  vibrational  and 
rotational  states  in  LTE  at  temperatures  T^.  and  ,  respectively,  is  also  calculated. 


The  first,  and  by  far  the  most  complicated  part,  is  the  distribution  over  the 
electronic  states.  The  MLTE  subroutine  is  best  described  in  terms  of  a  specific  example. 
Consider  the  following  system:  Gq  oxygen  atoms/cc  and  carbon  atoms/cc  at 

temperature  T,  =  T  ( ®K).  The  molecular  and  atomic  species  that  we  will  include  in  the 
calculation  are: 


1  G, 

2  Ol 

3  CO 

4  CO* 

5  Q 

6  O 

7  O* 

8  C 

9  C* 

10  e 

Note  1:  In  this  exan^le  and  throughout,  molecules  (in  arbitrary  order)  are  first,  then 
atoms  and  electrons  are  last. 

Note  2:  The  set  does  not  include  Cj .  If  there  are  any  stable  Cj  states,  they  are  not  listed 
in  Herzberg  and  will  be  ignored  here.  Also  CGj  is  not  included  since  the  code  is  limited 
to  diatomic  molecules. 


Next  let  us  define  a  descriptive  array  such  that  5,  *  is  the  number  of  atoms  of  type 

k  in  the  items  in  the  above  list.  It  is  convenient  to  extend  the  k  types  to  include  electrons. 
The  order  in  each  group  is  not  important,  however  the  code  is  hard  wired  to  have  the 
molecular  species  first  followed  by  the  atomic  species  and  then  the  electrons. 
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/ 


I, 

k=l  k=2 

C  O 


k=3 

e 


1  O2  0 

2  0 

3  CO  1 

4  CO^  1 

5  C2  2 

6  0  0 

10*  0 

8  C  1 

9  C*  1 

10  ^  0 


2 

2 

I 

1 

0 

1 

1 

0 

0 

0 


0 

-1 

0 

-1 

0 

0 

-1 

0 

-1 

1 


We  will  follow  the  procedure  (and  to  some  extent  the  notation)  of  Chapter  5  of 
“Atomic  Theory  of  Gas  Dynamics”  by  Bond,  Watson  and  Welch  [Ref.  10].  The  partition 
function  (their  equation  5-121)  can  be  written  as 


z  =  s  n  (1) 

t 


For  the  product,  /  runs  over  all  the  species  included  in  the  list,  and  the  sum  is  over  all  the 
sets  of  N,  consistent  with  the  following  constraints: 


10 

=  Cc 

/=! 

Conservation  of  carbon  atoms 

(2) 

10 

Z'VA,  =  Go 
/=! 

Conservation  of  oxygen  atoms 

(3) 

10 

XwA.,  =  o 

/=! 

Conservation  of  charge. 

(4) 

4 


The  partition  function  is  maximized  when  the  set  N,  has  the  “most  probable’ 


value,  N,.  That  is 


The  Z,{N,)  can  be  written 


=  n  Z.('V,)  =  r 


z,(N,)= 


(5) 


(6) 


The  Q,  are  further  defined  as 

Q,  =  v 


4  kM  kt 


3/2 


(7) 


where  the  Q,j.  are  the  sums  over  the  internal  states  of  the  species,  and  will  be  discussed  in 

subsection  Ha.  The  energy  levels  and  statistical  weights  of  the  internal  states  are  part  of 
the  input  files.  Note  that  for  electrons  that  sum  is  2. 


(8) 


By  using  Sterling’s  approximation  we  can  write 

10 

ln(r)  =  £  {-W,  [ln(Ar, )  - 1] + W,  •  ln{a )} . 

The  change  m  ln(r)  due  to  an  aifaitrarv  change  in  the  set  N,  is 

10 

5ln(r)=  £5^,  \n{QJN,), 

/=1 

however  we  need  to  subject  the  variations  &V,  to  the  constraints  of  eqs  2,3,  and  4,  Thus 


(9) 


/=! 

10 

L  /=! 

10 

£  N  ,S,^, 


/=l 


10 


0  -4 

£  =  0 

/=! 

10 

(10) 

0  ^ 

£  =  0 

/=! 

10 

(11) 

0  -> 

£  SN  ,S,  ,  =  0  . 

/=! 

(12) 

With  the  constraints  in  this  form  we  can  use  the  Lagrange  multiplier  technique  to  eliminate 
the  constraints.  That  is 
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(13) 


10  3  f  10  1 

«  ln(r)  =  X  [5N,  ■  \n{Q,  /  /V,)}  +  ^  A,  ^  SN,S„ 

'  =  1  t=1  I  /=!  J 

and  now  since  the  6N,  are  arbitrary  we  can  equate  the  coefficient  of  each  to  zero 
independently.  Thus  for  each  species  in  the  list  we  can  write 

ln{Q,  /  N  ,}+  Y,  A  *5,  ,  =  0  .  (14) 

*  =1 

It  is  convenient  to  define  n,  =  N,  /V  and  q,  =Q,  IV  then  eq.  14  can  be  written 


ln(^,  /  «,)  +  X  ^  =  0  .  (15) 

*  =  1 

To  eliminate  the  Lagrange  multipliers  we  note  that  there  are  three  rows  of  the  S  matrix 
that  have  a  single  entry  equal  to  1.  In  the  present  example,  these  rows  correspond  to  /  = 
6,  8,  and  10  i.e.  to  neutral  carbon,  oxygen,  and  the  electrons.  By  inspection  then 

A  I  =  -  ln(? g  /  /Jg ), 

A  2  =  -  ln(^6  /  ),  (16) 

A  3  =  -  ln(^,„  /  n,o  ). 

When  these  are  substituted  into  eq.  15  we  obtain 

ln(^,  /  n, )  -  ln(^g  /  n,  )S,  ,  -  ln(q,  /  n,  )S,  ,  -  ln(q,„  /  /!,„  )S,  ,  =  0  (17) 

or 

ln(^,  /«,)  =  ln{(^,  /«8r  'U,o  /«.o  )"'•’}  2  (18) 

«/  =  9/(«6  ’(«8  /^8f'(«,o  (19) 

Equation  19  can  be  used  to  solve  for  the  number  density  for  species  1, 2, ...,  9  in  terms  of 
two  unknowns  tig  and  /ig.  The  electron  density  /!,<,  can  be  calculated  from 

9 

^10  ~  ^1  9  (AO) 

/=i 

where  n^and  rtg  are  still  unknown. 

It  is  convenient  to  define 
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(21) 


=  R  «8  /  ^8 

V  2  =  ^  «  6  /  <7  6 

}'3  =  R  'Uo  /  ^10 

where  R  is  a  scaling  parameter  to  be  defined  at  the  time  of  numerical  solutions.  In  tenns 
of  these  new  variables  equation  19  becomes: 

3 

'h  =  ^/  n  (y*  /  R)  .  (22) 

*=i 

where  again  is  an  implicit  function  of  and  y2  ■ 

To  find  the  functional  dependence  of  y^  on  y,  and  we  write: 

9 

^10  ~ 

1=1 

=-E  «,n[3’./«]"s.., 

1=1  it=i 


Since 


=  -S  b./«]  s,., 

1=1  lr=l 

s  1 

Sj  3  is  either  -1  or  0,  [ J 3  /  7?  ]  Sj  3  =  - 5,-  3  and  we  can  write 

y^  /  R 


_  .=1  t=l 

”10 


therefore 


y^l  R 

3^3  =  R/Jio  / ^,0  =  [R  / ^,0]  — ^ 


or 


^>’3 = /^r‘-5,-.3  /^,o 

1=1  jfc=i 


}l/2 

■ 


(23) 


We  will  consider  two  cases:  constant  number  densities  and  constant  pressure. 
For  constant  number  density  we  can  rewrite  equations  2  and  3  as 
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(24) 


-Gc  =0 

/=1 
9 

y!  fh^i.z  -  ~  0- 

/=! 

When  equations  22  and  23  are  substituted  into  equations  24  we  obtain  two  equations  in 
two  unknowns. 


For  the  case  of  constant  pressure  we  need  to  insure  that 

to 

Ya't,  -  P^=0.  (25) 

One  more  equation  is  needed  to  close  the  system.  The  code  is  set  up  for  the  “no 
diffusion”  approximation.  In  this  approximation  we  can  write 

A  ® 

^  ^i,\  ~  ^co  ^  ^1,2  ~  0  *  (26) 

/=1  /=! 

where  is  the  constant  ratio  of  oxygen  to  carbon  nuclei. 

It  is  not  clear  that  either  equation  25  or  26  is  separately  consistent  with  the 
constraints  imposed  on  the  variation  procedure  used  to  find  the  most  probable  distribution 
(eq.  2, 3,  and  4).  One  can  imagine  an  iterative  scheme  in  which  constant  pressure  is 
achieved  by  first  solving  for  the  LTE  distribution  under  constant  number  density  and  then 
adjusting  the  number  densities  to  both  obtain  the  desired  pressure  (eq.  25)  and  number 
density  ratio  (eq.  26).  In  that  sense  eq.  25  and  26  together  may  satisfy  the  constraints  of 
eq.  2, 3,  and  4.  The  constant  pressure  version  of  the  code  was  implemented  replacing  eqs. 
24  by  25  and  26  in  the  numerical  scheme.  The  results  were  checked  against  the  constant 
number  density  version  and  found  to  be  in  remarkable  agreement. 

Equations  24  or  25  &  26  are  nonlinear.  On  the  other  hand,  the  Sj  j  are  integers 
and  in  fact  are  only  0,  + 1, 2 .  The  systems  of  equations  are  ideal  for  a  numerical  scheme 
such  as  the  BROYDN  subroutine  described  in  Chapter  9  of  “Numerical  Recipes  in 
Fortran”  by  Press,  Teukolsky,  Vetterling,  and  Flannery  [Ref.  13] .  For  the  special  case  of 
a  single  homonuclear  molecule  the  system  of  equations  reduces  to  a  single  equation  and 
brent’s  method  is  more  appropriate  then  BROYDEN’s. 
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Ila.  Sum  over  States 


The  first  point  to  be  addressed  is  the  zero  or  reference  point  for  bound  energies. 
Chemist  like  to  work  with  the  Enthalpy  and  to  have  the  ground  state  of  certain  molecules 
be  the  reference  point.  When  this  is  done  one  needs  to  include  the  heat  of  formation  in  a 
consistent  manner.  We  will  adopt  the  simpler  strategy  of  letting  the  ground  state  of  atoms 
be  the  reference  point  of  all  (bound  state)  energies.  With  this  reference  system  the  ground 
state  of  atomic  ions  has  an  energy  equal  to  the  ionization  potential,  while  the  lowest  state 
of  a  molecule  has  an  energy  equal  to  minus  the  dissociation  (into  the  ground  state  of  the 
atoms)  potential.  All  other  electronic  states  of  a  molecule  are  referenced  to  the  lowest 
state. 


We  look  first  at  the  sum  over  states  for  the  atoms  and  atomic  ions.  An  isolated 
atom  has  an  infinite  number  of  states  and  the  sum  over  them  diverges!  An  atom  in  a  gas 
or  plasma  is  not  isolated  and,  although  we  can  not  solve  the  lO"  lO'®  body  problem  to 
prove  it,  the  number  of  bound  states  is  finite  and  not  too  large.  For  atoms  in  a  plasma  the 
problem  has  been  addressed  by  Theimer  and  Kepple  [Ref.  14]  among  others.  Although 
similar  arguments  can  be  made  for  the  case  of  an  atom  in  a  gas,  we  will  adopt  a  much 
more  pragmatic  point  of  view  and  simply  truncate  the  set  of  bound  states.  From 
experience  we  know  that  the  results  are  not  sensitive  to  the  exact  number  of  bound  states 
that  are  retained.  Suffice  it  to  say  that  a  state  should  not  be  retained  if  its  mean  radius  is 
larger  than  that  of  the  unit  sphere.  That  is 


xl/3 


4;r  n 


(27) 


\  a  j 

where  is  the  density  (in  cm~^)  and  is  the  mean  radius  of  the  bound  state.  This  cutoff 
is  density  dependent  and  should  not  be  considered  an  upper  bound  to  the  cutoff  rather 
than  one  to  use  in  a  typical  calculation. 


The  sum  over  states  for  atoms  is 

^/  =  Z  Si  exp(-E,.  I  kT),  (28) 

f 

where  the  sum  is  over  the  set  of  bound  states  for  the  1  th  atom,  g.  is  the  statistical  weight, 
and  E.  =  0  for  its  ground  state.  For  atomic  ions  the  sum  over  states  is  similar 

=  Z Si  exp(-(/,  +E.)  I  kT)  =  exp(-/,  /  kT)'^  g.  exp(-E.  /  kT)  (30) 

*  I 

where  I,  is  the  ionization  potential  of  the  atom  of  the  1  th  species  and  again  E,  =  0  for  the 
ionic  ground  state. 
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Tlie  sum  over  states  for  molecules  is  similar  except  that  the  summation  is  over 
electronic,  vibrational  and  rotational  states. 


<?,  =exp(+D,/<T)  X  «,■ 

/,  V,  / 

In  the  above  equation  D,  is  the  dissociation  energy  from  the  ground  state  of  molecules  of 
the  1  th  species. 

The  basic  data  for  molecules  is  not  known  to  great  precision.  The  potential  curve 
for  each  electronic  state  is  characterized  by  only  a  few  parameters,  and  many  of  these  are 
only  valid  near  the  bottom  of  the  well.  The  energy,  in  wave  numbers,  is  given  by 
Hei^berg  [Ref.  1 1  equations  IV,  9;  IV,  10:  and  IV,  1 1]; 

T  =  r+G  +  F,  (32) 


G  =  o)^{v  + 1  /  2)  -  o)^x^{co  + 1  /  2)^  +  (o^y^iv  + 1  /  2)^+..., 


F  =  R/(/  +  1)-d/V  +  1)^+ .... 


(33) 

(34) 


e  e  e 


P  =D 
e 


8ty  jT  5  a  cP'g)  \ 


e  J 


In  the  above,  v  is  the  vibrational  quantum  number  and  J  is  the  rotational  quantum  number. 
In  spite  of  the  fact  that  the  parameters  were  fitted  to  the  bottom  of  the  potential  curve,  we 
retained  all  temis  such  that  the  energy  was  less  than  the  dissociation  energy  G  and  F 
greater  than  zero,  ■<  0.4  D, ,  and  1  -  DJ{J  + 1)  ^  0.1 .  These  cutoffs  in  V  and  J  are 
arbitrary  but  necessaiy.  Blind  use  of  equations  32  -  34  can,  and  often  does,  lead  to 
energies  that  are  negative  with  respect  to  the  ground  state  zero  point!  There  is  a  large 
body  of  literature  on  corrections  and  approximations  for  this  molecule  or  that  one.  We 
embarked  on  this  project  to  see  just  how  far  we  could  get  with  no  more  data  than  that 

found  in  [Ref.  1 1  and  12].  Thus  our  results  will  not  always  agree  with  those  from 
CHEMKIN. 
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III.  USER’S  GUIDE 


Subroutine  MLTE  calculates  the  LTE  distribution  over  the  electronic  states.  Once 
the  distribution  over  the  electronic  states  is  known  a  call  to  ALTEDIST  or  MLTEDIST 
will  calculate  the  distribution  over  the  excited  atomic  states  or  the  vibrational  and 
rotational  molecular  states. 

The  input  arguments  passed  to  subroutine  MLTE  are  the  electron,  vibrational  and 
rotational  temperatures  aU  in  degrees  Kelvin.  The  LTE  (electronic  states)  population 
distribution  is  returned  in  the  array  den.  The  populations  are  ordered  as  in  the  S  array 
(specified  in  the  input  file  “mlte.in”). 

The  atomic  and  molecular  data  are  read  by  MLTE  on  the  first  call  only.  In 
addition  to  the  formal  arguments  in  the  MLTE  call,  information  is  passed  to  and  from 
MLTE  through  the  following  common  blocks:  qp,  cntrl,  names,  and  logic.  The  use  of 
the  variables  in  these  common  blocks  is  illustrated  in  the  sanqile  MLTE  driver.  To  make  it 
easier  to  find  these  variables  they  are  printed  in  bold  face.  Many  of  these  variables  are  set 
via  the  namelist  “mltein”.  Default  vales  for  most  of  the  variables  in  the  namelist  are  set  in 
the  sample  driver. 

Some  important  variables: 
q  see  eq.  15  and  next  3  lines, 
en  population.. .see  eq,  22. 

r  numerical  scaling  factor. . .  see  eq.  2 1 . 
pO  (constant)  pressure, 

tk  electron  temperature  in  degree  Kelvin, 

fii  see  Fc0  in  eq.  26. 

amass  mass  of  the  atom  or  molecule. 

fg  used  to  calculate  a  “first  guess”  of  the  atomic  number  densities, 
dnl  total  (nuclear)  density  of  first  atomic  species. 

xlow/xhigh  For  the  case  of  one  homonuclear  molecular  species,  xlow/xhigh  are  the 
min/max  number  density  for  the  neutral  atom, 
ztol  Tolerance  in  the  solution  of  the  equation  for  the  neutral  atomic  number 

density  in  the  above  case. 

eps  Convergence  tolerance  for  the  case  of  more  than  one  unknown, 
imol  Number  of  molecular  species.  Specified  in  “mlte.in”. 
iatom  Number  of  atomic  species.  “ 

id  imol  +  iatom 

idp  id  +  1 

maxloop  maximum  number  of  iterations. 
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index  For  internal  use  only! 

is  Tlie  “S”  array.  Specified  in  “mite. in”, 

ns  Number  of  independent  variables, 

nsp  ns  +  1 

jw  Number  of  vibrational  levels  included  in  the  calculations, 

jvr  Number  of  rotational  levels 

cp  Constant  pressure? 

latom  Calculate  the  distribution  over  the  excited  atomic  states? 

“  vibrational  and  rotational  molecular  states? 

Ibug  Print  debug  information? 

Iplot  Non-standard  Plots? 

Inp  Use  xiow/xhigh? 


The  arguments  to  ALTEDIST  are  the  electronic  state  distribution  (output  from 
MLTE)  and  the  electron  temperature.  The  arguments  to  MLTEDIST  are  the  electronic 
state  distribution,  the  electron,  vibrational,  and  rotational  temperatures.  At  present  the 
outputs  are  files  “mlte.atom”,  and  “mlte.mol”.  The  calls  to  ALTEDIST  and  MLTEDIST 
can  be  implemented  by  setting  latom  and  Imol  =  .trae.  before  the  call  to  MLTE. 

If  the  code  returns  a  “Singular  Jacobian”  error,  you  need  a  better  first  guess. 

If  the  code  returns  an  error  for  too  many  iterations  in  Broyden,  you  can  continue,  but  if  it 
the  error  messages  persist,  you  may  need  a  better  first  guess  for  the  number  density,  i.e. 
den(i).  In  the  sample  driver,  the  input  parameter  dnl  sets  the  overall  density  and  can  be 
set  in  the  namelist. 

The  gradients  in  the  partition  function  are  so  steep  that  unless  you  have  a  fair  “first 
guess  the  code  can  go  into  a  “chaotic  mode”.  The  Broyden  package  is  quite 
sophisticated  and  we  have  not  been  able  to  modify  it  to  avoid  the  fall  into  chaos.  From 
eq.  21  and  22  we  note  that  the  y’s  determine  all  the  n’s  but  only  the  n’s  corresponding  to 
the  neutr^  atomic  species  determine  the  y’s.  Thus  it  is  clear  that  it  is  best  to  start  at  some 
intermediate  temperature  where  our  intuition  tells  us  that  most  of  the  populations  will  be 
neutral  atoms  and  thus  we  can  make  a  very  good  first  guess  of  the  y’s.  Once  we  have 
obtained  a  converged  solution  at  any  temperature,  we  can  then  “adiabatically”  (the 
converged  solution  at  the  last  temperature  is  a  good  first  guess  for  the  next)  change  the 
temperature  up  or  down  as  needed.  The  same  technique  may  be  necessary  for  very  small 
pO,  small  or  large  fh  (ratio  of  atomic  species). 
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Ilia.  EXAMPLES 


Tlie  following  four  examples  all  use  the  sample  driver  listed  on  pages  29  -  33.  For 
these  examples  the  desired  output  are  plots  of  the  populations  Vs  temperature.  The  files 
♦.corn  are  UNIX  scripts.  Their  purpose  is  to  show  explicitly  the  input  and  data  files. 


INPUT  DATA 
casel.com 

cp  /home/mydisk/lte/casel.m  inlte.iii 

/hoine/mydisk/Ite/casel.nlst  mlte.nlst 

/hoine/mydisk/data/molq/tnolqo2.in  molel 

/bonie/mydisk/data/inolq/molqo2p.in  inole2 

yhome/mydisk/data/molq/molqco.in  mole3 

/home/mydisfc/data/molq/molqcop.in  inole4 

/bome/myclislc/data/mol^molqc2.in  moleS 

/honie/mydiskydata/inolq/atqo.in  atoml 

/hoine/niydisk/data/molq/atqop.in  atom2 

cp  /home/mydisk/data/molq/atqc.in  atom3 

cp  /home/mydisk/data/molq/atqcp.in  atom4 

/home/mydiskAte/mlte.x 


easel. nist 

$mltein 
cp=.tnie., 
p0=i.3332d3, 
fn(l)=1.0, 
fn(2)=1.0, 
itag='C,0  1.0’, 
fg=0.001, 
diil=1.0el2, 
r=1.0dl6, 
scale(l)=1.0, 
scale(2)=  100.0, 
scale(3)=1.0, 
eps=0.001, 
maxloop=80, 
npp=4, 
t0=0.690, 
delte  =  -0.000431, 
ntp=850, 
latom=.false., 
lmol=.faIse., 

Ibugl=.false., 

Ibug2=.false., 

Ibug3=.false., 
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ltrace=.false., 

jvv=65, 

jvr=300, 

tvv=:-1.0, 

trr=-1.0, 

dslip(  1  )='$$$$$$$$$$$$$$$$$$$$$$', 

dshp(2)='$$$$$$$$$$$$$$$$$$$$$$'^ 

dshp(3)=’$$$$"$$"$$$$$"$$"$', 

dshp(4)='$$"$"$’'$$$"$"$"$', 

dshp(5 

dshp(6)='$$$$$$$$$$$"$$"$$"$', 

dshp(7)='$"$"$"$$$"$"$"$$’, 

dshp(8)='$"$$"$$"$$""$$"$', 

dshp(9)=’$$"$$$"$$$$"$$$"$$', 

nlb(l)=28, 

nlb(2)=28, 

nlb(3)=28, 

nlb(4)=28, 

nlb(5)=28, 

nlb(6)=28, 

nlb(7)=28, 

nlb(8)=28, 

nlb(9)=28, 

$end 


cp . constant  pressure? 

pO . if  cp  is  true,  then  the  pressure  =  pO 

ft* . ratio  of  atoms  (nuclei)  of  ty|>e  n  to  typ>e  1 

fg . used  for  first  guess 

r . scale  factor 

scale . more  scale  factors 

eps . convergence  test 

maxloop...max  iterations  for  convergence 

laton . print  out  atomic  population  dist  ? 

Itttol . print  out  molecular  population  dist  ? 

Ibug . print  out  debug  info  7 

jvv . number  of  vibrational  levels 

jvr . number  of  rotational  levels 

Ivv . vibrational  temp  (  ®K)  if  negative  tvv  =  tk 

trr . rotational  temp  (  ®K)  if  negative  trr  =  tk 
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easel. in 


5 


1  02 
2  02  + 

3  CO 

4  C0+ 

5  C2 

6  O 

7  0+ 

8  C 

9  C+ 
10  e 


4 

S 

k=l  2  3 

C  O  e 
0  2  0 
0  2-1 
110 
11-1 
2  0  0 
0  10 
0  1-1 
10  0 
10-1 
0  0  1 


In  this  example  it  was  necessary  to  start  at  an  intermediate  temperature  (8000 
and  step  down  to  lower  temperatures  and  then  make  a  second  run  to  work  up  to  higher 
temperatures.  The  namelist  for  the  second  ran  was  the  same  as  for  the  first  except  for  the 
sign  of  delte. 

In  the  following  figures  the  As  are  from  CHEMKIN  and  indicate  the  agreement 
between  MLTE  and  CHEMKIN.  Note:  The  CHEMKJN  data  base  is  only  valid  for 
temperatures  below  6000  ^K. 
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10.0ITORR)  C.O  1.0  mute  H.T.  P0*  100.0(TORR) 


TFMPFRATURe  (1000  K)  TEMPERATURE  (1000  KJ 


casel.coni 


cp  /homc/mydisk/ltc/case2.in  mltc.in 

cp  /home/mydisk/lte/casc2.nlst  mlte.nlst 

cp  /home/m ydisk/data/molq/molqo2.in  mold 

cp  /homc/niydisk/data/molq/molqo2p.in  mole2 

cp  /home/niydisk/dala/molq/molqno.in  niolc3 

cp  /home/mydisk/data/molq/molqnop.in  mole4 

cp  /home/mydisk/data/molq/molqn2.in  mole5 

cp  /home/m ydisk/data/molq/molqn2p.in  niole6 

cp  /home/niydisk/data/molq/atqo,in  atom  1 

cp  /home/mydisk/data/molq/atqop.in  atom2 

cp  /home/mydisk/data/molq/atqn.in  atom3 

cp  /home/mydisk/data/molq/atqnp.in  atom4 

/home/mydisk/lte/mlte  .x 


casel.nist 


$mltein 
cp=.true., 
p0=1.3332d4, 
fn(l)=1.0, 
fh(2)=1.00, 
itag='  N,0  1.00', 
fg=0.001. 
dnl=1.0el2, 
i^l.0dl6, 
scale(l)=1.0, 
scale(2)=100.0, 
scale(3)=1.0, 
eps=0.001, 
maxloop=80, 
npp=4, 

10=0,069, 

delte  = -0.000650, 

ntp=850, 

latom=.false., 

lmoI=.false., 

lbugl=.false., 

Ibug2=,false., 

Ibug3=.false., 

Itrace=.false., 

lplotl=.false., 

lplot2=,  false., 

x2low=3.80d9, 

dlx2=1.05, 

tpl=+8.0, 

dtp=0.04, 

spl=5.0, 

dsp=0.02, 
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jvv=65, 

jvp=300, 

tvv=-l.O, 

trr=-1.0, 

dshp(  1  )=’$$$$$$$$$$$$$$$$$$$$$$', 

dshp(2)=’$$$$$$$$$$$$$$$$$$$$$$', 

dshp(3)='$$$$"$$"$$$$$"$$"$'. 

d.shp(4)=$$"$"$"$$$"$"$”$', 

dshp(5 

dshp(6)='$$$$$$$$$$$"$$"$$"$', 

dshp(7)=$”$"$"$$$'T$"$$’, 

dshp(8)=’$"$$''$$"$$"’'$$"$’, 

dshp(9)='$$"$$$"$$$$"$$$"$$', 

nlb(l)=28, 

nlb(2)=28, 

nlb(3)=28, 

nlb(4)=28, 

nlb(5)=28, 

nlb(6)=28, 

nlb(7)=28, 

nlb(8)=28, 

nlb(9)=28, 

Send 


cp . constant  pressure? 

pO . if  cp  is  true,  then  the  pressure  =  pO 

fii . ratio  of  atoms  (nuclei)  of  type  n  to  type  1 

fg . used  for  first  guess 

r . scale  factor 

scale . mote  scale  factors 

eps . convergence  test 

maxloop...max  iterations  for  convergence 
laton.....ptint  out  atomic  population  dist  7 

Imol . print  out  molecular  population  dist  ? 

Ibug . print  out  debug  info  ? 

jw . .number  of  vibrational  levels 

jvr. . number  of  rotational  levels 

tw . vibrational  temp  ( ‘’K)  if  negative  tw  =  tk 

trr . rotational  temp  (  "K)  if  negative  trr  =  tk 
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casc2.m 
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1 
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10.0{TORRI  N,0  1.00  MLTE  H.T.  P0*  100.0(TORR)  N.O  1.00 


ON 


10.0tTORR)  N.O  1.00  mute  M.F.  P0*  100.0ITORR)  N.O  1.00 


TDtPERATURE  (100B  K)  •  ■  TEMPERATURE:  (1000  K1 


case3-coni 


cp  /home/mydiskAte/case3.in  mlte.in 

cp  /home/m ydisk/lte/case3.nlst  mlte.nlst 

cp  /home/mydisk/data/molq/molqn2.in  mold 

cp  /home/mydisk/data/molq/molqn2p.in  mole2 

cp  /home/m ydisk/data/molq/atqn. in  atoml 

cp  /home/mydisk/data/molq/atqnp.in  atom2 

/home/mydisk/lte/mlte.x 


case3.nlst 


$mltein 
cp=.true., 
p0=L0000e3, 
fn(l)=1.0, 
fia[(2)=L00, 
itag=’N2,N  \ 
fg=0.00001, 
dnl=1.0el6, 
r=1.0dl6, 
scale(l)=1.0, 
scale(2)=100.0, 
scale(3)=L0, 
eps=0.001, 
maxloop=80, 
npp=4, 
t0=L207, 
delte  =  -0.000202, 
ntp=501, 
xlow=1.0e-36, 
xlugh=1.0e+16, 
latom=,false., 
lmol=.false., 
lbugl=.true., 

Ibug2=.false., 

Ibug3=.false., 

ltrace=:.tnie., 

lplotl=.fasle., 

lplot2=.false., 

tpl=-9.0, 

dtp=0.02, 

jvv=35, 

jvp=150, 

tvV=:-L0, 

tlT=-l.O, 

dshp(  1  )=*$$$$$$$$$$$$$$$$$$$$$$•, 
dshp(2)=’$$$$$$$$$$$$$$$$$$$$$$', 
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dshp(3)='$$$$''$$"$$$$$"$$"$', 

dshp(4)='$$"$"$"$$r$"$"$’, 

dshp(3 

d.slip(6)='$$$$$$$$5;$$"$$"$$"$', 

d.shp(7)='$"$"$"$$$"$”$"$$', 

d.slip(8)='$"$$"$$"$$""$$"$’, 

dshp(9)='$$"$$$"$$$$"$$$"$$', 

nlb(l)=28, 

nlb(2)=28, 

nlb(3)=28, 

nlb(4)=28, 

nlb(5)=28, 

nlb(6)=28, 

nlb(7)=28, 

nlb(8)=28, 

nlb(9)=28, 

$end 


cp . constant  pressure? 

pO . if  cp  is  true,  then  the  pressure  =  pO 

fa . ratio  of  atoms  (nuclei)  of  type  n  to  type  1 

fg . used  for  first  guess 

r . scale  factor 

scale . more  scale  factors 

eps . convergence  test 

maxloop...max  iterations  for  convergence 

laton . print  out  atomic  population  dist  ? 

Imol . print  out  molecular  population  dist  7 

Ibug . print  out  debug  info  ? 

jvv . number  of  vibrational  levels 

jvr . number  of  rotational  levels 

tvv . vibrational  temp  (  ®K)  if  negative  tvv  =  tk 

trr . rotational  temp  (  °K)  if  negative  trr  =  tk 


case3.in 
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75.0ITORR) 


TEMPERATURE  [1000  K1 


casc4.com 


cp  /honie/mydi.sk/ltc/casc4.in  nilte.in 

cp  /home/niydisk/lte/case4.nlst  mlte.nlst 

cp  /home/mydisk/data/niolq/molqn2.in  mold 

cp  /home/m ydisk/dat a/m olq/molqn2p. in  molc2 

cp  /home/mydisk/dala/molq/atqn.in  atom  1 

cp  /home/m ydisk/dala/molq/atqnp.in  atom2 

cp  /home/m ydisk/data/molq/atqar.in  atom  3 

cp  /home/mydisk/data/molq/atqarp.in  atom4 

/home/mydisk/lte/mlte.x 


case4.nlst 


$mltein 
cp=.true., 
pO=L3333e3, 
fn(l)=1.0, 
fn(2)=3.0, 
itag='  N2  N  At 
fg=0.00001, 
dnl=1.0e2, 
r=1.0dl6, 
scale(l)=l,0, 
scale(2)=I  000.0, 
scale(3)=1.0, 
eps=0.001, 
maxloop=80, 
npp=4, 
t0=1.2000, 
delte  =  «0.0030, 
ntp=350, 
xlow=1.0e-16, 
xhigh=1.0e+16, 
iatom=.false., 
lmol=  false,, 
lbugl=.true., 
Ibug2^.false., 
Ibug3=.false., 
ltrace=.true., 
lplotl=  fasle., 
lplot2=.false., 
tpl=-9.0, 
dtp=0.02, 
jvv=35, 
jvi^I50, 
tvv=-1.0, 
trr=-1.0, 
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dshp(l)=’$$$$$$$$$$$$$$$$$$$$$$', 

dshp(2)='$$$$$$$$$$$$$$$$$$$$$$', 

dshp(3)='$$$$"$r$$$$$"$$"$', 

dshp(4)=’$$"$"$"$$$’'$"$"$', 

dshp(5)='$"$"$"$"$"$"$"$’, 

dshp(6)=’$$$$$$$$$$$"$$"$$"$', 

dshp(7)='$"$"$"$$$"$"$"$$', 

dshp(8)=’$"$$"$$"$$""$$"$', 

dshp(9)='$$"$$$"$$$$"$$$"$$', 

nlb(l)=28, 

nlb(2)=28, 

nlb(3)=28, 

nlb(4)=28, 

nlb(5)=28, 

nlb(6)=28, 

nlb(7>=28, 

nlb(8)=28, 

iilb(9)=28, 

Send 


q) . .constant  pressure? 

pO . .if  cp  is  tnie,  then  the  pressure  =  pO 

fii . ratio  of  atoms  (nuclei)  of  type  n  to  type  1 

fg . used  for  first  guess 

r. . scale  factor 

scale . more  scale  factors 

eps . convergence  test 

maxloop...max  iterations  for  convergence 

laton . print  out  atomic  population  dist  ? 

Imol . print  out  molecular  population  dist  ? 

Ibug . print  out  debug  info  ? 

jvv . number  of  vibrational  levels 

jvr . number  of  rotational  levels 

tw . vibrational  temp  (  ®K)  if  negative  tw  =  tl 

trr . rotational  temp  (  ®K)  if  negative  trr  =  tk 
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la.BiTORRl  N2  N  Ar  HLTE  M.F.  P0*  100.0(TORR) 


TCTPERATU^  tieae  K)  .  TDIPERATURE  flBBB  K1 


IV.  SAMPLE  DRIVER 


cdeck  tmlte  < - «< 

PROGRAM  tmlte 
c 

c  driver  for  subroutine  mite... molecular  Ite  distribution 
c 

implicit  real*8  (a-h,  o-z) 
c 

c  from  MLTE5.MTR  September  2, 1997 
c 

integer  va,  vd,  vm,  vam,  vamp,  ve,  vap,  vdvm2,  w,  vr 
PARAMEnrER(va=10,  vap=va+l,  vd=20,  vm=10,  vam=va+vm, 

.  vamp=vam+l,  ve=30,  vdvm2=vd*vm*vm,  vv  =  50,  vr=100) 
logical  cp,  latom,  Imol,  Ibugl,  Ibug2,  Ibug3,  Iplotl,  lpIot2, 

.  Inp 

character*6  sname 

dimension  index(va),  is(vam,  vap),  sname(vamp) 
common  /cntrl/  eps,  imol,  iatom,  id,  idp,  maxioop,  index, 

.  is,  ns,  nsp,jvv,jvr 
common  /names/  sname 
save  /cntrl/,  /names/ 

common  /logic/  cp,  latom,  Imol,  Ibugl,  lbug2,  lbug3,  Iplotl, 

.  lplot2,  Inp 
save  /logic/ 

dimension  q(vamp),  en(vamp),  scale(va),  fn(va),  amass(vamp) 
common  /qp/  q,  en,  r,  pO,  tk,  scale,  fn,  amass,  fg,  dnl, 

.  xlow,  xhigh,  ztol,  tpl,  dtp,  spl,  dsp,  x21ow,  dlx2 
save  /qp/ 
c 

common  /bugs/  Itrace 
c 

logical  Itrace 
c 
c 

parameter  (NPPTS=4,  NTPTS=1000,  NSM=20) 
c 

character*51  titl 
character*9  itag 
character*28  dshp(NSM) 
character*  16  agdshn 
c 

real*4  temp,  dens 
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dimension  temp(NTPTS),  dens(NTPTS,  NSM),  nlb(NSM) 
data  titl/’MLTE  M.F.  P0=  (TORR)  $'/ 
dimension  den(vamp),  sden(v;unp,  500),  sfk(500) 
c 

namelist  /mltein/  cp,  pO,  dnl ,  fn,  fg,  r,  scale,  eps,  latom, 

.  Imol,  Ibugl ,  lbug2,  lbug3,  jvv,  jvr,  dnl ,  tvv,  trr, 

.  maxloop,  ntp,  npp,  xlow,  xhigh,  ztol,  Iplotl ,  lplot2,  tpl,  dtp, 
.  tO,  delte,  pO,  itag,  dsbp,  nib,  Itrace,  x21ow,  dlx2,  spl,  dsp 
c 
c 
c 

call  opngks 
c 

c  default  values  for  namelist  variables 
c 

cp  =  .true, 
iatom  =  .false. 

Imol  =  .false. 

Ibugl  =  .false. 

Ibugl  =  .false. 

Ibug3  =  .false, 
xlow  =  l.Oe-16 
xhigh  =  l.Oe+16 
ztol  =  l.Oe-06 
Itrace  =  .false. 

Iplotl  =  .false. 

Iplotl  =  .false, 
r  =  1.0el6 
eps  =  l.Oe-4 
fg  =  0.90 
tvv  =  -1.0 
trr  =  -1.0 
jvv  =  10 
Jvr  =  50 
ntp  =  NTPTS 
npp  =  NPTPS 
to  =  0.01 
delte  =  0.002 
p0=  1.0d4 
maxloop  =  70 
x21ow  =  l.Oe-16 
dlx2  =  3.0 
spl  =  -5.0 
dsp  =  0.02 
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3 


do  3  i  =  I,  va 
iiide,\(i)  =  0 
scale{i)  =  1.0 
fn(i)  =  0.0 

open(unit=4,file='mlte.nlst',status=’old') 
open(unit=6,  file='mlte.out',  status='unknown') 
open(unit=7,file='mlte.tl',status='unknown') 
open(unit=8  ,file='mlte.t2'  ,status='unknown' ) 
open(unit=9,f ile=’nilte.b  r,status='unknown' ) 
open(unit=22,  file- mite,  tbl',  status- unknown') 
c 

c  read  namelist 
read(4,  mltein) 
write(6,  mltein) 
c 

if  ((.not.  Iplotl)  .and.  (.not.  lplot2))  then 
c 

call  agseti('DASH/SELECTOR.’,9) 
call  agseti('DASH/LENGTJ.',  28) 
callagsetf(’DASH/DOLLAR-QUOTE.',0.004) 
do801i=l,9 

call  agsetc(agdshn(i),  dshp(i)) 

801  continue 
c 

call  agsetc('LABEL/NAME.',  'L') 
call  agsetiCLINE/NUMBER.’,  100) 
call  agsetc('LINE/TEXT.’,'MOLE  FRACTIONS’) 
c 

call  agsetc(’LABEL/NAME.’,  •B')  ' 
call  agseti(’LINE/NUMBER.',-100) 
call  agsetc('LINE/TEXT.';$') 
call  agseti(’LINE/ISrUMBER.’,-120) 
call  agsetc('LINEA^XT.’, TEMPERATURE  (1000  K)$') 
c 

endif 

c 

c 

c 

do  200  np  =  1,  npp 
Inp  =  .tme. 
pO  =  10.0*  pO 
c  write(6, 913)  pO 
913  fomiat( '  pO  =  ’,  lpel2.3/ 
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.  '  temp  02  02+  CO  CO+', 

C2'/'  O  0+  C  C+') 

titl(32:40)  =  itag 
pOl  =  pO/ 1. 3332(13 
write(titl(  16:22),  99)  pOt 
99  fonnat(n.l) 
write(6,  92)  pO 

92  fonnatC  Pressure  =',  lpel2. 3) 
c 
c 
c 

te  =  tO 

do  100  nt  =  1,  ntp 
te  =  te  +  delte 
tk  =  te*  1 1600.0 
temp(nt)  =  tk/l  000.0 
tvv  =  tv 

if  (tvv  .le.  0.0)  tvv  =  tk 
trr  =  tr 

if  (trr  .le.  0.0)  trr  =  tk 

c 

c 

call  inlte(  tk,  tvv,  trr,  den) 
c 

if  (nt  .eq.  1)  call  nilte(  tk,  tvv,  trr,  den) 
if  (nt  .eq.  1)  call  nilte(  tk,  tvv,  trr,  den) 
c 
c 

write(6, 90)  te,  tk,  (den(i),  i  =  1,  id) 

90  format(  ‘  te  =  ’,  f8.4,  ’  (eV) '  tk  =’,  f8.0, '  (K)’/ 
.  Cden=Mp4el5.4)) 
c 
c 

stk(nt)  =  tk 
do  23  i  =  1 ,  idp 
23  sden(i,  nt)  =  den(i) 
c 

sumden  =  O.OdO 
do  13  i  =  1,  id 

13  sumden  =  sumden  +  den(i) 
do  12  i  =  1,  id 

c  pden(i)  =  den(i)/sumden 
12  dens(nt,  i)  =  den(i)/sumden 
c 
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c  write(8, 912)  tk,  (pden(i),  i  =  1,  NSM) 
c  912  fonTiat(/  Ip5el5.4/  (Ip5el5.4)) 
c 
c 

100  continue 
c 

write(7,  9122)  temp 

9122  fomiatC  temp  =’/  (Ip8el0.2)) 
do  9123  ip  =  1,  id 

9123  write(7, 9124)  (dens(it,  ip),  it  =  1,  ntp) 

9124  format(lp8el0.2) 

call  ezmxy(temp,  dens(l,  1),  NTPTS,  id,  ntp,  titl) 
c 
c 

write(22, 925)  pO 

925  formate///  Pressure  = lpel2.3) 
do310j=  l,ntp 
c 

xntot  =  0.0 
do320ii=  l,idp 
320  xntot  =  xntot  +  sden(ii,  j) 
c 

do  330  ii  =  1,  idp 
330  sden(ii,  j)  =  sden(ii,  j)/xntot 
c 

write(22, 923)  stk(j),  xntot,  (sden(i,  j),  i  =  1,  idp) 
310  continue 

923  fonnat(fl0.2,  elO.3,  8fl0.7) 
c 
c 
c 

Inp  =  .false. 

200  continue 
c 

call  clsgks 
c 

stop  333 
c 
c 
c 

END 


33 


oncocofooofooooo 

Doooooooooo 

c^<T><T>a^o^a^ooo»H 

C^CNlCMCSICMC^tOlOlOU) 


rH<NfHO<T> 
CMrHO^OO  10t0l0U?<NJ 

M<V|COCOOOCJCM<MiO<MC'4 


VDU5U>U>Vr)lOV£>VD 

OrHrHrHrHrHrHrHr-f 

Doooooooooocoooco 

lomioioioioiotn 


TT 

r'  r-  <r>  o  CO  fH 

(DlOlOCsIO^^lOiHCh 

i^OOlOVDOOCOCn<Tk 

csjcsjrorooocviroco 


iHt-li-IHr-fTHr-lrH 


lO  r-l  ^  U>  00 

^  ^  <n  o  o  o  <M  ID 
J^r-tD^lD  V>00^<Ti00 
Oir-liHi-fHCJiHt-ICNirKCO 
HOOOOOOOOOO 

oooooooooo 


^  ^  <n  lo  00 
43oo<t^ooiocm^ 
(X<j^<T»ooi^r^oooor- 

WOOOOOOOO 

oooooooo 


oo^m<n  r^^cic'jr* 

0)  <7\  D4  CO  ^  CM  CO  CM  00  a\ 
CQrMOOCnrMI^IOOOC^OOOO 
oocD^CDoor^r~rMoor^ 


<U  CM  CM  MO  00  U>  ID 
CQ«HH^OrOrMlOr^ 
OkO^^OOr-fr-tO 
^'DCOcnCMVOCMCM 

THrHiHrHf-MrHrMtH 


i 


> 


i0 

+» 

<d 

Q 

U 

<d 


3 

Cl 

o 

H 

I 


ooooooo  oo 

OOOOOOOf-lOO 

ooooooor-oo 

OOOOOOOCMOO 

NOOOOOOOtHOO 

oooooooooo 

+++++++ I ++ 


rf  ^ 

inocooooooo 

OIOIOM>0000 

OCMCMrMOOOO 

NOOOOOOOO 

^oooooooo 

oooooooo 

+  +  4-  I  I  +  +  + 


o 


o 

I 


r- 

o 

00  o 

CD 

to 

CM  f-H 

o 

o 

O  O 

to 

00 

•  « 

• 

• 

o  o 

o 

CM 

+  1 

1  + 

+  + 

00 

CM 

o 

+ 


I’ 


lO  rH  to  ^  ^  O 

oiocy\oor-CDioo 
f-f^Or'rMlDlOO 
OOOOrHr'OO 


UO 


OOOOOOOO 

I  +  I  1111  + 


r- 

00 

o 

o 

to 

I 


o  o  o  o 
OOtOrHOr-^OOCMCT^O 
COrMHCMCOtDIOa^rOO 


'OK*HrH00lOt£>  rj* 

^  >:ootDtDnooooo 

o  CM<n^tDu>^r-cM 


I  ^  O  O  O  <Tk  O  O 

r-IrMrMrMrMrHrHrM 


rH  to  to  to  CM  to  r* 

r'CO^COtOCMrMlOtOtO 

>^i+OOOHOO<7ktO<JkrM 

tO'M’r^otooooocMr- 

OOVD'crtDCTkr-OOrHOOVD 


o 


c 

•H 


•  ^  o  o  to  “M* 

^  OCMCMOiHtOCOtO^ 


<UOtO<M««-fCM<MiHtO<Tk 
Eh  t~fC0<Tkr^CMtDC\CO 
r-  ^cocoocMr^cM 
VO  00  CO  O  ^  O  CO 

rH  CM  CO  rr  Tf 


:s 


Cl 

tr 

1  d 

0 

(3^  tn 

tn 

a 

+ 

04 

1 

+ 

+ 

o 

CO 

CO 

cu 

to 

04  04 

04 

to 

tr^ 

<l) 

00 

H 

rH 

00 

tH 

00 

00  rH 

00 

rH 

2 

§ 

00 

E  a 

:z; 

X 

X> 

< 

o 

d 

C 

4) 

Q 

r- 


rj* 

00 

o 

to 

to 


D 

+ 

CO 


o 


O  rHrHOrrCM^CM 
00  ^  to  <Tk  CM  CM 


O>C0  00rHr-00CMTr 
'O^frHr-rHrHOO 
rHf^CMrHrHtOOO 

CM  rH  rH  rH  rH  rH  rH  rH 


00  O  ^  00  TT 

CM  O  to  rH  CM  •  •  • 

jy . to  ^  OO 

EHOVOtOOOr^OOCM 
00  CM  CM  oo  O  O  <Tk 
CD  00  rH  CM  to  to  to 
CO  to  rH  ^  CD  CO  VjD 
q  TT  to  to  CD 


o 

V 

tr 


go  S 

6oz 


+  +  I  I 

COCUCOQCOOiCOQ 

rHOOCOOOCOrHrHTH 


34 


molqcop . in 

C0+  28.0  112917.0  -89462 


in 


• 

00 

• 

• 

r- 

r- 

CM 

r4 

<Ti 

cr» 

CM 

CM 

O 

o 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

<J\ 

r- 

CO 

CO 

CO 

CO 

o\ 

O 

Q 

CD 

CO 

CO 

CO 

CO 

rr 

o 

CO 

CO 

in 

CM 

CO 

00 

00 

00 

m 

in 

o 

OO 

00 

00 

00 

00 

00 

00 

00 

00 

00 

Q 

<0 

00 

r- 

r- 

<Ti 

a\ 

o\ 

C' 

r- 

o\ 

o 

O 

o 

o 

Q 

<y\ 

<Tt 

r- 

<Ti 

r- 

<y\ 

r- 

r- 

<Tt 

o 

rH 

rH 

t— 1 

rH 

rH 

rH 

rH 

rH 

rH 

<j\ 

<y\ 

Ok 

r- 

r- 

r- 

00 

CTi 

CM 

o 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CM 

00 

in 

m 

00 

O 

OO 

r- 

o 

CO 

00 

CO 

CO 

<y\ 

CM 

rH 

CO 

in 

00 

00 

d) 

00 

r4 

r- 

CM 

CO 

H 

rH 

rH 

<D 

rH 

r- 

|-«- 

v>| 

ID 

t— 1 

tH 

CM 

CO 

r4 

CO 

u 

o 

CM 

CM 

CM 

CM 

CM 

CM 

CO 

rH 

rH 

in 

rH 

rH 

o 

o 

CM 

I  t-l  H  H  rH  H 


H  t-|  H 


^  CM  ID 
a>  ^  CM 
OiOO  CTk  O  ^ 
H  H  fH  00  CM 
id  O  O  O  O 
•  •  •  • 
o  o  o  o 


00 

tH  rH  ^  r-  00  CO 

COOCT^^lOCkU^H^ 

j^r‘cor-^u>r-co\0oo 

tH^^r^r^r^r^^^HCMCM 

HOOOOOOOOOOO 

ooooooooooo 


H  00 
CO  00 

JCj  00  00  ^  ^  o 

(XH  H  CM  rH  CM  in 

H  o  O  O  O  O  O 

<d  «  •  •  •  •  « 

o  o  o  o  o  o 


fl)  CM  ^ 

05  o 
r*  00  <Ti  in 
in  00 


^  in  00 

C4  \0  ^  CO  c\  o\  r^ooiD 

<l)oo^r^oo<J^iDcooo^r^<7> 
ooc^inoor-f^HOCMCMC^^ 
a^^iorr^co^cnoocTko 


in  00 
0)  CM  ^  ^  r-l  00  CTk 

cQ  00  Tj»  T-i  o 
ov  o  o  tH  in 

i-f  rM  CM  CM  r~l  ,H 


O  O  O  O 
O  O  O  O 
N  O  O  O  O 

>  o  o  o  o 
o  o  o  o 

o  o  o  o 

+  +  4-  + 


r^  CM  CO  0> 

CMO^IDOOCMCMOOOOO 

oHoor^ooooinoo 

MOOOOOOOOOOOO 

^ooooooootnoo 

ooooooooooo 

I  l  +  l  l+H’+l+M- 


o  o  o  o  o 
o  o  in  o  o  o 

N  O  O  O  O  O 
>00^000 
o  o  o  o  o  o 

o  o  o  o  o  o 

+  4-  I  +  +  -h 


«H  00  o 
O  CO  00  o 
O  T-l  CM  O 

O  O  00  o 
•  «  •  « 
o  o  o  o 
14-4-4- 


tn 

>iVo  CO  cn  00  00 

>CMOO<J^OOCM<J^OOOO 
CMOlOrMf-tr-OOCO 
OrHtn^^OOOOO 

ooooooooo 

OOOOOOOOCM 

I  4-  I  4-  4-  4-  I  14- 


*1  o  in  o  o 

:  O  O  O  CO  O 
^  O  00  O  c~l  00 

o  o  m  o  o 

•  ••••• 

o  o  o  o  o  o 

I  I  I  1  I  I 


I  ^  r* 

:  CO  00  CM 

fH  in  cn  00 

♦  •  •  • 

m  00  r-  00 

t4  rM  CM  00 


ri 

r- 

00 


r*  rH 

!  ^  CM  CM  rM  ^  0>  m 

cMr-cMooi^^ooin^ 

OOOOrMr-IOC^COOO^ 

^00M*<MCM00f-lr-l00 

iHr-ltHr-CrMiHr-ltHCM 


^  i 

^  O  CM 

o\ 

rH 

<7^ 

in 

r-  J 

S  rH  O 

rH 

<J> 

CM 

00 

•  • 

• 

•  • 

COV) 

00  ^ 

H 

C\ 

00 

00 

1  * 

rH  rH 

CM  rH 

rH 

^  CO  00 
CM  O  tH 

I  •  •  •  • 

Xt  CM  ^  'M* 
rM  CO  00  ^ 
CM  in  r-  rH 

CM  rH  t-l  fH 


00  00 
r-'M‘cytooinocoa»r- 
incoOOOOCMCMCM^r-t 

ooooocooM»c\cMi^int-i 

lncooor^ooo^ln^^ooa\ 

oo^r-tnincoini^or-ir- 

CMrHi-ltHtHf-lt-l  CMCM 


lo 

CO 


O  O 

o  r-  00 


r-  in 


'o  >r-ooo%<x>r-fH 

m  o  o  t-l  <y>  o 

CM  CM  ^  00  <T»  O 
1-1  CM  rH  CM  CM  CM 


00  r- 

o  •  .  . 

M)  •  00  CO  CM 

B  o  00  tH 
r-  00  o 
o  in  00 

CM  ^  CO 


00 

CM 


00  in  r-  CO  CO 
O  CO  00  -M*  CO  o  "M* 


00 

00 


<DOOOCnCMCMOOr^OCDOOr-l 

E-c  OtHr-mcoovoromin 

CMCOCMtHCMOO\rHOOOO 
o  o>  CO  00  o^  CM  r-  o^  in  00 
inmcococor^oooochcyk 


oo  < 

CM  <u 
M 


a 

•H 


in 

in 

CM 

rH 

<j\ 

• 

• 

« 

« 

• 

rH 

r- 

00 

OO 

CO 

CO 

co 

rH 

o 

CO 

rr 

00 

CO 

rH 

m 

in 

CM 

o\ 

CM 

CM 

in 

CO 

-i-  4- 

co  cu  cn  Q 

O  CM  CM  CM  CM 

I 

Jz;  X  «<  pQ  o 


4-4-11  4- 

cocnoicocoatQQOLicoat 


gtHOOOOOOtHrHtHOOCOOOOO 
O  CM  S  -  - 

u  ^  ^  X  <  0^  cQ  id  id  ^oawo 


CM 


04 

CM 

& 

H  . 

2  ^ 

fl  :z; 


d  tn  d 

4-  4-4-  + 

CO  ou  CO  CO  cu  CO 

<I)  CM  CM  CM  rr  CM  CM 

§ 

5s  X  rtj  CQ  id  Q  O 


35 


orrooor-VDvo 
QVDcgogcMior-r- 
CM  m  in  m  lo  cTk 
CM  -M*  TJ*  xr  00  o  o 
lo  u>  VO  v£>  CO  r' 


rocovovocovovovo 

OiHi-fVOVOCDVOVOVD 

OcTkOkcocococoron 

COOO<T»a^Ch<T»<7\<Ti 


CMCMCMCMCMCMCOOO 

r^t^r-r'r'r-tooo 

O<T^<TiCr»<T>cn<y»Tff0CM 

Ooooooovoiovo 


00  rr  00  »-i 

O  CO  VO  CM  H  VO 

fl)invoiHvovovovo 

Mr-IO’^OOOO 
fH  fH  tH  1-!  H  H  fH 


lO  to 
00  rH  CM 

J^I^O^CMOHOOCO 
CUtH  r-f  1-1  CO  CM  i-f  CM 
HOOOOOOO 

rt . 

o  o  o  o  o  o  o 


VO  tn  VO  CO 

«)Tf«>CMOCMVOCM 
(Qocninoooooo 
o  T^  o  o  cn  o\ 

r-l  i-t  fH  CM  (M  iH  iH 


O  CM  r-l 

>^CM'^tr>oor^ooco 
cooor^oocnoooocTv 
VOCMHCMCMCMCMH 
. . 


USO^CM^CMCMCMCMtO 
OiOO  00  CT>  Ch  <T^  Tf 
H00f-li-trHi-tr-lf-tCM 
(di-IOOOOOOO 

. . 

•  ooooooo 
o 


a)r-<Tk<^r^r-cOiHvo 

CQo^vDcor^mvovooo 

<7>(OU>COCO<OC040 


CM  lO 

ID  IT)  r'  ID 

r-  lo  VO  r-  iH  rr  co 

a)OHCMi-fOOCMO’^rH 

MCMCMCMinTfinVOHTH 


CO  O  iH  VO  VO 

tM  CM  c\  HO 

jCini^OOCOVD^CMVOO 
CUH  HHHCMHHCMCM 
HOOOOOOOOO 
ld*****«*«a 

ooooooooo 


VO  r- 

fl)VO^CO  VOCT^OCO 

rOlOVOOlO  OHCOO 

^cMonvonoor-r- 

asj*  ^  ^  o  <y>  •  •  # 

. .  H  H 

H  H  H  O  O  O 


OOOOOOO 
OOOOOOO 
NOOOOOOO 
^OOOOCMOO 
O  O  O  O  CM  O  O 

OOOOOOO 
+  +  +  +  I  +  + 


OOOOOOOO 

NOOOOOOOO 

^OOOOOOOO 

OOOOOOOO 

OOOOOOOO 

OOOOOOOO 

++++++++ 


r^ooioooooo 

NCMOOIOOOOOO 

^HOOOOOOOO 

OOOOOOOOO 

OOOOOOOOO 

OOOOOOOOO 

I  4-  +  +  +  +  +  +  + 


iO  lO 
H  VD 

H  ^  r*  to 

O  O  H  1^ 

•  •  •  • 

OOOOOOO 
+  I  +  +  +  +  + 


>1  00 

>  io 

OOOOOOOCM 

OOOOOOOO 

+++++++ I 


O 

<n 

00 

O 

o 

o 

o 

CM 

in 

H 

o 

o 

• 

in 

• 

• 

• 

•  O  O  o 

O  I  I  I 


a\ 

m 

a\ 

• 

00 

• 

2 

1 

0 

6 

7 

5 

o  r- 

<j\ 

in  VO 

00 

<n 

VO . 

•  o 

c\ 

H  o 

•  • 

CO 

c\ 

XCMtO^OOCOVOCsl 

00 

i  O  H  CO 

in  CM 

CO 

t  •  H  H  H  H  H 

CM 

m  9  m 

H  CM 

•  • 

CM 

VO 

o 

in 

XT  VO  00 

VO  o 

m 

H 

CM 

1 

H  H 

H  CM 

1 

o 

o 

CO 

• 

00  H 

o 

• 

CM  O  O  O  O  O 

o 

o 

O  00  00 

o\ 

CO 

CO 

•M*  •  •  •  •  « 

•  •  « 

•  • 

•  • 

o 

•  CO  O  CO  Tf  CO  00  H 

?^^o>incoin^ 

CM 

^  VO  <7)  H  H  00  00 

r-  o 

o  r-  CO 

0\  CM 

r-  <j\ 

r-  CM  r-  00  CM  CM 

CM  VO 

cn  00  o  CO  CO  CO  CO 

r- 

CO  H  H  H  H  H 

H  CM  H 

CM  CM 

CM  CM 

CM 

o 

o 

CM 

• 

0 

7 

6 

0 

7 

8 

• 

o  o  o  o  o  o 

o 

o 

•  •  • 

«  • 

•  • 

o 

o  lO  CM  VO  ^  00  O 
<y  VO  Tf  CM  00  CM  O 

EH  CT^  0\  H  o  VO  00 

CO  lO  CM  n  o  H 
^  ^  m  in  VO  VO 


OOOOOOOH 
I  CT^  rr  00  CM  ^  to  r' 

I  H  CM  00  r-  in  ^  rr 

CM  <7>  rH  Oh  rH  CO 

lO  in  VO  VO  VO  r- 


<T\  VO 

•  00  o  00  VO  in 

0\  ovcTkor*  Hvor-o 

lo  H  •  •  •  . . 

CM  ^HCM^CMOCMOO^O 
H  rHHHHCMHHHH 

I 


a\  r-  CT> 

H  in  CM 


^ococM^oo^cnr-r* 

ooooco<y^ina\otncM 

tn’^^c^oor^t^cncn 


►  ocoinr^or-cooH 
0)  •HChincricy>cy\vocr> 
Enocni-Hovocor-cMO 
r-coco^tn<T>invo 
H  CO  CO  CO  r-  r- 


4  +  4  4 

CO  CU  CU  CO  CO  Q 

OCMCMCMCMCMCMCM 


isxcmoaMPu 


O  4  4  I  4 

C  COCOOiOCOCOOCu 

^  0)  H  CO  CO  CO  00  rH  rH  rH 

oSS 


b^tT>b>:3:3t7':3dd 

^  I  +1  + I +  + 

O  COQCOCOQCOCOCOCO 

tr 

^  gCOHHHOOCOCOCOH 

o  CM  i 

6o:z:x*«u:)  oict:ri:oQQVM 


36 


molqo2p . in 

02+  32.0  96786.0  -41259 


^  ^  O 

o  VO  ^  lo 

Q  <M  CM  CM 
CM  CM  CM  r- 
m  IT)  IT)  VO 


00 

00  O  VO 
<u  VO  tH  <Tk  <T\ 
M  rM  00  O  r* 

1-!  00  CM 


VO  to  VO  VO 

jC  r-  00  o 

CuO\  lO  CM 
H  r-l  rM  T-l  CM 
«J  O  O  O  O 

O  O  O  O 


VO  0\ 
00  VO  CM 

0)  (Ti  O  VO  00 
(0  VO  iH  O  CM 


O  O  O  O 
O  O  O  O 
O  O  O  O 
N  O  O  O  O 
>  O  O  O  O 

o  o  o  o 
+  -f  +  + 


o  o  o  o 
o  o  o  o 
o  o  o  o 
o  o  o  o 
o  o  o  o 


o  o  o  o 
+  +  -f  + 


to  cTt  o 
CM  00  to  o 


VO  O  00 


olor- 
r-vocMr- 
^  • 
to  00  VO 
o  00  ov  <n 

CTk  O  00  T-( 
rH  r-t  tH 


O  00  O 
o  •  •  • 

0)  •  ^  <J>  CM 

Em  o  VO  VO  m 
VO  to 

CM  O  ov 
00  -M* 


ty>  rJ  0 

I 

0^  Oi  Pm  CO 
<U  CM  -M*  CM  'M* 


S  X  <d  rij  XJ 


37 


Atomic  Data 
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VI.  SPIXIFIC  HFAT 


The  input  to  CHEMKIN  is  the  Enthalpy  and  its  derivative  witli  respect  to 
temperature,  .  We  start  by  calculating  the  excitation  energy: 


E.g,  exp(-E,.  I  kT) 

fJ,  =  - - •  (35) 

^g,  exp(-£,  /kT) 

I 

In  Eq.  35  the  summation  is  over  bound  states  for  atoms  and  over  electronic,  vibrational 
and  rotational  states  for  molecules.  Ionization  and  dissociation  potentials  are  not  needed 
since  they  cancel  in  the  numerator  and  denominator.  At  this  point  we  also  define: 

X^XexpC-E,./^?’) 

W,  =  ^ - .  (36) 

'^g.Qxp{-E.  /kT) 

I 

From  Eq.  35  we  can  calculate 

^  =  (W,-u;)/kT\  (37) 

By  definition: 

H,  =  Hf+^kT  +  U,,  (38) 

thus 

^  ^  f  +  (W^/  -  ^' )  /  (kTf  ].  (39) 

The  results  of  these  calculations  are  presented  in  Figures  25-48.  These  results 
are  presented  first  as  the  cumulative  function  as  each  electronic  state  is  added.  In  the  next 
figure  the  MLTE  results  (solid  line)  are  compared  to  those  of  [Ref.  7]  (dashed  line). 
Recall  that  the  data  of  [Ref.  7]  is  only  valid  up  to  6000  °K.  For  all  atoms  and  all 
molecules  other  tlian  and  the  MLTE  results  are  a  little  higher,  but  very  close  (up  to 
6000  °K)  to  those  of  [Ref.  7].  Also  note  that  for  all  molecules  other  than  Cj,  the 
cumulative  is  a  monatonic  nondecreasing  function  of  electronic  states.  That  is,  each 
electronic  states  has  a  positive  contribution  to  the  total  for  the  molecule.  For  Cj,  the 
<3 

addition  of  the  “a  first  excited  electronic  state  to  the  ground  state  contribution, 

significantly  decreases  .  Each  additional  electronic  state  contributes  positively  to  . 
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TruDro /( -rrior  t  *  nnf^  i/>  i 


Cum-dU/ dT 


TEMPERATURE  (1000  K)  TEMPERATURE  (1000  K) 


With  the  technique  described  in  Appendix  A  the  heat  capacity  data  can  not  only  be 
put  into  the  convenient  and  compact  form  of  a  fourth  order  polynomial  in  T,  it  can  also  be 
made  to  be  continuous  to  the  data  of  [Ref.  7]  at  1000  °K.  In  this  form  the  heat  capacity 
data  is  shown  in  Table  I  and  is  valid  over  a  temperature  range  from  1,000  °K  to  20,000 
°K. 


TABLE  I 

4 

Coefficients  for  the  expansion:  C  =  X  X.T 

P  ;  _  n  * 


02 

02+ 

NO 

NO+ 

N2 

N2+ 

CO 

CO+ 

C2 

o 

0+ 

N 

N+ 

C 

c+ 


1,000°K  ^  T  <  20,000°K 


^0 

3.84564E+00 
4.30658E+00 
3.81893E+00 
4.95483E+00 
4.92123E+00 
3.04607E+00 
4.9973  lE+00 
4.07416E+00 
3.85173E+00 
2.51193E+00 
2.69824E+00 
2.61047E+00 
2.51903E+00 
2.38164E+00 
2.49994E+00 


X^ 

3,84552E-04 

-3.54956E-04 

2.95232E-04 

-1.31580E-03 

-1.26664E-03 

1.05458E-03 

-1.33305E-03 

-2.32754E-04 

6.16230E-04 

-5.80432E-06 

-2.68845E-04 

-1.85537E-04 

-2.69791E-05 

1.38052E-04 

6.28564E-06 


X, 

-3.49093E-08 
1.65847E-07 
-2.66967E-08 
3.03675E-07 
2.88429E-07 
-1.24073E-07 
3.45703E-07 
1.48826E-07 
-1.2403  lE-07 
1.16068E-08 
7.56701E-08 
8.19666E-08 
1.43068E-08 
-2.03215E-08 
-5.11023E-09 


X3 

2.39363E-13 

-1.48660E-11 

5.76425E-13 

-1.43951E-11 

-1.23135E-U 

5.81310E-12 

-2.18446E-11 

-1.25583E-11 

7.67680E-12 

-1.16638E-12 

-5.18263E-12 

-7.08549E-12 

-1.12102E-12 

1.44489E-12 

8.67188E-13 


^4 

2.41539E-17 

3.70888E-16 

-2.66467E-19 

1.15619E-16 

5.11080E-17 

-1,05991E-16 

3.87341E-16 

2.89677E-16 

-1.56928E-16 

3.51696E-17 

1.08000E-16 

1.82018E-16 

2.66681E-17 

-2.78334E-17 

-2.55265E-17 


From  equation  39  we  see  that  the  Enthalpy  can  be  calculated  from  the  heat 
capacity  data  by  integration.  Thus 


H(j)  =  \l^Cpdr  =  + 


H{T)  = 


(40) 

(41) 
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In  [Ref.  7]  the  reference  temperature,  T„,  is  taken  to  be  298.15  °K.  We  will  use  this 
reference  and  their  data  for  up  to  T,  =  1,000°K  to  evaluate 


H\T,)  = 


In  the  notation  of  [Ref.  7]: 


H\T^)  =  R 


i  =  4 
b,+  1 
1  =  0 


fl,r, 


i+I 


1  +  1 


(42) 


(43) 


The  coefficients  b^  and  a,,  are  tabulated  in  [Ref.  7]. 


CHEMKIN  is  set  up  to  use  the  data  in  the  same  format  as  that  of  [Ref.  7].  Thus 
we  need  to  calculate  y,,  such  that  we  can  have 


H  iT)  =  R 


i  =  4  X.  T*' 
3'.+  I  -T— - 
/  =  0  '  +  ^ 


(44) 


It  is  clear  that 


yi=bi  + 


and  are  tabulated  in  Table  H. 


(45) 


TABLE  II 


3’i 

y^ 

^>2 

-1.295342E+03 

0^, 

1.394564E+05 

^2 

1.795155E+05 

^2 

1.708964E+04 

CO 

-1.512556E+04 

CO^ 

1.486719E+05 

NO 

9.693373E+03 

NO^ 

1.17351  lE+05 

<^2 

9.902393E+04 

0 

2.9243 12E+04 

0^ 

1.878475E+05 

N 

5.606134E+04 

2.256235E+05 

C 

8.550795E+04 

2.168825E+05 
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We  can  now  use  the  data  from  Tables  I  and  II  to  replace  the  CHEMKIN  data  for  T 
>  1,000  °K,  and  thus  extend  the  validity  of  the  CHEMKIN  results  up  to  20,000  °K. 
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APPENDIX  A 


LEAST  SQUARE  EIT 
WITH  A  CONSTRAINT 

We  wish  to  fit  a  set  of  data  {x,  ,  /  }  to  a  polynomial  in  x.  It  is  known  that,  if  the 

range  of  x  is  large,  the  usual  algorithm  results  in  a  nearly  suigular  set  of  linear  equations. 
To  avoid  this  problem  we  introduce  a  set  of  basis  polynomials  O,. .  Ideally  this  set  of 
polynomials  are  orthogonal  over  the  range  of  jt,..  That  is 

b 


where  a<x.<b.  For  simplicity  we  will  use  a  set  of  polynomials: 


<E>, 

=  "1.1 

=  ^2.1 

+ 

^2.2 

a: 

% 

=  "3.. 

+ 

^3,2^ 

+ 

^4 

=  "4,, 

+ 

^4,2 

X 

+ 

«4.3^' 

+  "4.4^' 

^5 

=  ^5.I 

+ 

^5,2 

X 

+ 

+  aj4A-'  +a^^x* 

(Al) 


and  select  the  coefficients  for  approximate  orthogonalization.  The  set  of  coefficients  is 
listed  in  Table  I  and  plots  of  Oj  -  Oj  are  shown  in  figures  AI  -  A4.  It  follows  that  if 

F(jc)  =XC’y"’,then 

1=1  1-1 

C3  =  ^5  <35  5 

C4  =  ^3  4  + 

C3  —  B^  ^3  3  +  B^  ^4  3  +  B^  ^3  3  (A2) 

Cj  =  B^  a^  2  ^4  ^4.2  ^3  ^3,2  ^2  ^2,2 

“  ^5  ^5,1  ^4  ^4,1  ^3  ^3.1  ^2  ^2.1  ^1,1  • 


In  the  (unconstrained)  least  square  algoritlim  one  would  minimize  the  residual 


"|2 


/=! 


(A3) 
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If,  in  addition  to  minimizing  the  above  residual,  one  needs  to  impose  the  constraint 


=  =  (A4) 


then,  by  the  technique  of  Lagrange  multipliers,  one  can  minimize  the  new  function: 


2H[  ]. 


(A5) 


This  minimization  is  accomplished  by  setting  the  derivatives  of  R  with  respect  to  the 
coefficients  R,.  to  zero.  Thus 


+  2A®,{a:-)  =  0, 


i=t  *  * 

for  j  =  1,  2,  ...,5.  Next  the  equation  for  j  =  1  can  be  solved  for  A: 


A  =  -! - 


5 

s 


On  substitution  of  equation  A7  into  equation  A6,  one  obtains 


1=1  k 


-E/* 


0,(x‘) 


(A6) 


(A7) 


=  0,(A8) 


But  O,  is  a  constant,  therefore  equation  A8  simplifies  to 


(A9) 


for  j  =  2,  3, 4,  and  5.  These  four  equations  plus  equation  4  constitute  a  linear  set  of  five 
equations  in  the  five  unknown  B,..  If  the  Os  are  near  enough  to  orthogonalality,  the  linear 
equations  will  be  well  behaved. 
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The  linear  set  can  be  written  VB  =  W ,  where  V. ,  =  =  f\ 


Kj  -  and 

*  *' 

forj  =  2,..,5. 


(AlO) 

(All) 


TABLE  AI 

Coefficients  of  C>,. 


^M  = 
«2.1  = 
«2.2  = 
«3.I  = 
^3,2  “ 
«3.3  = 
«4.I  = 
^4,2  ~ 
^4.3  = 
«4.4  = 
"^5,1  = 
^5,2  “ 
«5.3  = 
"5.4  = 
"5,5  = 


+  1.0 
-1.0 
+1.0e4 
+1.0 
-2.67e-4 
+1.30e-8 
-  1.0 
+4.13e-4 
-4.68e-8 
+1.56e-12 
+1.0 
-4.88e-4 
+7.76e-8 
-4.96e-12 
+1.10e-16 
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